      program mxsim25h0                            !  H0 first 100 data
      implicit real*8 (a-h,l,o-z)
      real*8 lb,lb2
      logical max
      parameter(np=3,mp=np+1,np2=7,mp2=np2+1,
     &        nplay=80,nsim=4*nplay,ndec=100,ntype=5,
     &        ftol=1.0d-10,neps=4,nloop=50,nobs=4*nloop)
c                                          
      external famoeb1,famoeb2,gamoeb1,gamoeb2,fcn1,gcn1
      character*24 fdate
      dimension p(mp,np),x(np),y(mp)
      dimension p2(mp2,np2),x2(np2),y2(mp2)
      dimension gm(nsim),u5(nsim),u10(nsim)
      dimension pa(ndec,4),pb(ndec,4),ma(4),mb(4)
      dimension d1(ndec),d2(ndec),d3(ndec)
      dimension is1(ndec)
      dimension lb(np),ub(np),c(np),vm(np),xopt(np),xp(np),fstar(neps)
      dimension nacp(np)
      dimension lb2(np2),ub2(np2),c2(np2),vm2(np2),xopt2(np2),
     &        xp2(np2),nacp2(np2)
      dimension xzll21(nplay+1,nobs),rho(nplay),srho(nplay)
      dimension jlow(nplay+1),jhigh(nplay+1),in(nplay)
      dimension xlow(nplay+1),xhigh(nplay+1)
      dimension axz11(nplay),sxz11(nplay),axz21(nplay),sxz21(nplay)
      dimension bll(nobs),sxz2(nplay),s2(nplay)

c
      common /ch/is1
      common /mats/pa,pb,ma,mb,d1,d2,d3
c
      open(unit=10,file='eut75a.data')
      open(unit=11,file='HOparm.data')
c1
      do 2 n=1,ndec
         read(11,*)pa(n,1),pa(n,2),pa(n,3),pa(n,4),
     &       pb(n,1),pb(n,2),pb(n,3),pb(n,4)
         d1(n) = pa(n,2) - pb(n,2)
         d2(n) = pa(n,3) - pb(n,3)
         d3(n) = pa(n,4) - pb(n,4)
    2 continue     
c
      do 3 i=1,nsim
         read(10,*)gm(i),u5(i),u10(i)
    3 continue
c
      ma(1) = 0
      ma(2) = 10
      ma(3) = 20
      ma(4) = 30
      mb(1) = 0
      mb(2) = 10
      mb(3) = 20
      mb(4) = 30
c
      call tools
c-----------------------------------------------------------------
      write(*,'('' ***** mxsim25h0d ******'',20x,a24,
     &//"    HO Experiment, first 100"/)')fdate()
      write(*,'(10x,"Data simulated from EUT model."/)')
c              
      do 6 id=1,nplay+1
        jlow(id) = 0
        jhigh(id) = 0
        if (id.le.nplay) then
           axz11(id) = 0.d0
           axz21(id) = 0.d0
           sxz11(id) = 0.d0
           sxz2(id) = 0.d0
           sxz21(id) = 0.d0
           in(id) = 0
        endif
        do 5 nl=1,nloop
          xzll21(id,nl) = 0.d0
    5   continue
    6 continue
c    
      dseed = 90001.d0
      do 200 nl=1,nloop
       do 190 it=1,4
        mt = it + 4*(nl-1)
        do 120 i0=1,nplay
          id = it + 4*(i0-1)
          gam=gm(id)
          u1=u5(id)
          u2=u10(id)
          call genx(dseed,gam,u1,u2)
c
c     Estimate the Toobox model on the is1 simulated data
c Set parameters for simulated annealing
      max = .true.
      eps = 1.d0
      rt = 0.5
      iseed1 = 1947
      iseed2 = 4012
      ns = 20
      nt = 25
      maxevl = nt*(nt+5)*ns*np2/10 + 2
c      maxevl = 1
      iprint = 0
c
      do 20 i = 1, np2
         lb2(i) =  -5.d0
         ub2(i) =  5.d0
         c2(i) = 2.0
   20 continue
         x2(1) = dlog(1.d0/0.4d0 - 1.d0)
         x2(2) = dlog(1.d0/0.6d0 - 1.d0)
         x2(3) = dlog(1.d0/0.6d0 - 1.d0)
         x2(4) = dlog(1.d0/0.8d0 - 1.d0)
         x2(5) = dlog(0.2d0/0.05d0 - 1.d0)
         x2(6) = dlog(0.15d0/0.05d0 - 1.d0)
         x2(7) = dlog(0.10d0/0.05d0 - 1.d0)
c
      T = 8.d0
      do 30, i = 1, np2
         vm2(i) = 0.5d0*(ub2(i)-lb2(i))
   30 continue
c
      call sa(np2,x2,max,rt,eps,ns,nt,neps,maxevl,lb2,ub2,c2,iprint,
     &        iseed1,iseed2,t,vm2,xopt2,fopt,nacc,nfcnev,nobds,ier,
     &        fstar,xp2,nacp2,gcn1)
c
         do 35 k=1,np2
            x2(k)=xopt2(k)
   35    continue
c
c initialize starting vertices of the simplex. 
         delta=0.5d0
         call jnit(x2,p2,y2,delta,gamoeb1)
c 
c call the simplex routine. 
         ndim=np2 
         iter=3000
c         iter=1
         call amoeba(p2,y2,mp2,np2,ndim,ftol,gamoeb1,iter) 
c 
c         write(*,'(a,i5)')'iterations: ',iter 
c compute average values of vertices. 
         call avg(p2,x2,np2,mp2) 
         zll11 = -gamoeb1(x2)   
c
c compute LL21
         zll21 = -gamoeb2(x2)
c
c   Estimate the EUT model on the simulated is1 data.
c Set parameters for simulated annealing
        max = .true.
        eps = 1.d0
        rt = 0.5
        ns = 20
        nt = 25
        maxevl = nt*(nt+5)*ns*np/10 + 2
c        maxevl = 1
        iprint = 0
        do 50 k = 1, np
           lb(k) =  -5.d0
           ub(k) =  5.d0
           c(k) = 2.0
   50   continue
c
         gam = 0.01d0 + 0.98d0*gam/120.d0
         uh = 0.01d0 + 0.97d0*u1
         du = 0.01d0 + 0.97d0*(u2 - u1)
         x(1) = dlog(1.d0/gam - 1.d0)
         x(2) = dlog(1.d0/uh - 1.d0)
         x(3) = dlog((0.995d0-uh)/du - 1.d0)
c
        T = 8.d0
        do 60, i = 1, np
           vm(i) = 0.5d0*(ub(i)-lb(i))
   60   continue
c
        call sa(np,x,max,rt,eps,ns,nt,neps,maxevl,lb,ub,c,iprint,
     &       iseed1,iseed2,t,vm,xopt,fopt,nacc,nfcnev,nobds,ier,
     &       fstar,xp,nacp,fcn1)
c
1000  format(/,' simulated annealing summary',/,
     &       /,' initial temp: ', g8.2, '   rt: ',g8.2, '   eps: ',g8.2,
     &       /,' ns: ',i3, '   nt: ',i2, '   neps: ',i2,' maxevl: ',i10)
1001  format(/,' optimal function value: ',g20.13
     &       /,' number of function evaluations:     ',i10,
     &       /,' number of accepted evaluations:     ',i10,
     &       /,' number of out of bound evaluations: ',i10,
     &       /,' final temp: ', g14.6,'  ier: ', i3)
c
         do 65 k=1,np
            x(k)=xopt(k)
   65    continue
c
c initialize starting vertices of the simplex. 
         call init(x,p,y,famoeb1)
c 
c call the simplex routine. 
         ndim=np 
         iter=3000
c         iter=1
         call amoeba(p,y,mp,np,ndim,ftol,famoeb1,iter) 
c compute average values of vertices. 
         call avg(p,x,np,mp) 
c         write(*,'(i3,3g13.5)')nl,zll11,xll11,zxll11(id,nl)
c compute LL21
         xll21 = -famoeb2(x)
c		 
         xz11 = xll11 - zll11
         xz21 = xll21 - zll21
         xzll21(id,mt) = xz21
         if (xz11 .lt. -3.89) then
           in(id) = in(id) + 1
           axz11(id) = axz11(id) + xz11
           axz21(id) = axz21(id) + xz21
           sxz11(id) = sxz11(id) + xz11**2
           sxz2(id) = sxz2(id) + xz21**2
           sxz21(id) = sxz21(id) + xz11*xz21
         endif
c
  120    continue
c        
  190  continue
  200 continue
c      stop
c
      do 225 i=1,nplay
         xobs = float(in(i))
           axz11(i) = axz11(i)/xobs
           axz21(i) = axz21(i)/xobs
           sxz11(i) = sxz11(i)/xobs - axz11(i)**2
           sxz2(i) = sxz2(i)/xobs - axz21(i)**2
           s2(i) = dsqrt(sxz2(i)/(xobs-1.d0))
           sxz21(i) = sxz21(i)/xobs - axz11(i)*axz21(i)
           rho(i) = sxz21(i)/sxz11(i)
           if (xobs .gt. 2.d0) then
             s21 = (sxz2(i)/sxz11(i) - rho(i)**2)/(xobs-2.d0)
           else
             s21 = 0.d0
           endif
           srho(i) = dsqrt(s21)
c
  225 continue
      write(*,'(/"  in  ",80i5)')(in(i),i=1,nplay)
      write(*,'("axz11 ",80g13.5)')(axz11(i),i=1,nplay)
      write(*,'("axz21 ",80g13.5)')(axz21(i),i=1,nplay)
      write(*,'("sxz2 ",80g13.5)')(s2(i),i=1,nplay)
      write(*,'(/" rho ",80g13.5)')(rho(i),i=1,nplay)
      write(*,'(" srho ",80g13.5)')(srho(i),i=1,nplay)
c      
      end
c***********************************************************************
      double precision function famoeb1(x)
      implicit real*8 (a-h,o-z)
      parameter (np=3,ndec=100)
      dimension x(np),px(ndec,2)
      dimension is1(ndec)
      common /ch/is1
c
      gam=1.2d2/(1.d0+dexp(x(1)))
      u1=1.d0/(1.d0+dexp(x(2)))
      u2=u1 + (1.d0-u1)/(1.d0+dexp(x(3)))
c
      call feut(gam,u1,u2,px)
c
      t=1.d0
      do 10 k=1,75
         j = is1(k)
         t = t*px(k,j)
   10 continue
      famoeb1 = -dlog(t)
      return                   
      end
c***********************************************************************
      double precision function famoeb2(x)
      implicit real*8 (a-h,o-z)
      parameter (np=3,ndec=100)
      dimension x(np),px(ndec,2)
      dimension is1(ndec)
      common /ch/is1
c
      gam=1.2d2/(1.d0+dexp(x(1)))
      u1=1.d0/(1.d0+dexp(x(2)))
      u2=u1 + (1.d0-u1)/(1.d0+dexp(x(3)))
c
      call feut(gam,u1,u2,px)
c
      t=1.d0
      do 10 k=76,ndec
         j = is1(k)
         t = t*px(k,j)
   10 continue
      if (t .gt. 1.d-300) then
        famoeb2 = -dlog(t)
      else
        famoeb2 = 7.d2
      endif
      return                   
      end
c***********************************************************************
      subroutine feut(gam,u1,u2,px)
      implicit real*8 (a-h,o-z)
      parameter (ndec=100)
      dimension px(ndec,2)
      dimension pa(ndec,4),pb(ndec,4),ma(4),mb(4)
      dimension d1(ndec),d2(ndec),d3(ndec)
      common /mats/pa,pb,ma,mb,d1,d2,d3
c
      do 20 n=1,ndec
         yn = d1(n)*u1 + d2(n)*u2 + d3(n)
         pn = 1.d0/(1.d0 + dexp(-gam*yn))
         px(n,1) = pn
         px(n,2) = 1.d0 - pn
c         write(*,'(i3,2g12.4)')n,yn,pn
   20 continue
c      stop
c
      return
      end
c***********************************************************************
      subroutine init(xint,p,y,famoeb)
      implicit real*8 (a-h,o-z) 
      parameter (np=3,mp=np+1)
      dimension q(np,np),p(mp,np),x(np),y(mp) 
      dimension xint(np),xx(np),xxx(np)
      external famoeb 
c 
      do 1 k=1,np
         xxx(k)=xint(k)
    1 continue
c
c      write(*,'(//''Initial Parameter Values:'')')
c      do 14 k=1,np
c         write(*,'(3x,g12.6,3x,g12.6)')xint(k),xxx(k)
c   14 continue
c
c      y0=-famoeb(xxx)
c      write(*,'(/''initial log-likelihood = '',g14.7/)')y0
c      stop
c
      p(1,1)=1.d0
      p(2,1)=-1.d0
      if (np.eq.1) goto 51
      do 5 k=2,np
         p(1,k)=0.d0
         p(2,k)=0.d0
    5 continue
c
      do 50 it=2,np
         f0=float(it)
         f1=(dsqrt(f0*f0-1))/f0
         f2=-1.d0/f0
c
         do 11 j=1,it
            do 10 k=1,np
               q(j,k)=p(j,k)
   10       continue
   11    continue
         do 15 k=1,np
            p(1,k)=0.d0
   15    continue
         p(1,it)=1.d0
         do 21 j=2,it+1
            p(j,it)=f2
            do 20 k=1,it-1
               p(j,k)=f1*q(j-1,k)
   20       continue
   21    continue
   50 continue
   51 continue
c
c Initial step size
      delta=0.05d0
c      write(*,'(''  stepsize = '',g12.3//)')delta
c 
      do 80 j=1,mp
         do 79 k=1,np
            p(j,k) = xxx(k) + delta*p(j,k)
   79    continue
   80 continue
c 
c Initialize y. 
      do 92 i=1,mp 
         do 91 j=1,np 
            x(j)=p(i,j) 
   91       continue 
         y(i)=famoeb(x) 
c     write(*,'(i2,9(1x,g12.4),g14.7)')i,(x(j),j=1,np),y(i)
   92    continue 
      return 
      end 
c***********************************************************************
      subroutine tools
      implicit real*8 (a-h,o-z)
      parameter (ndec=100)
      dimension pa(ndec,4),pb(ndec,4),ma(4),mb(4)
      dimension d1(ndec),d2(ndec),d3(ndec)
      dimension ymax(ndec),ymin(ndec),ymp(ndec),ypm(ndec)
      common /mats/pa,pb,ma,mb,d1,d2,d3
      common /ys/ymax,ymin,ypm,ymp
c
      do 50 n=1,ndec
         maxa = 0
         mina = 15
         pa0 = 0.d0
         pma = 0.d0
         za0 = 0.d0
         maxb = 0
         minb = 15
         pb0 = 0.d0
         pmb = 0.d0
         zb0 = 0.d0
         do 10 j=1,4
            if ((pa(n,j).gt.0.d0).and.(ma(j).ge.maxa)) maxa = ma(j)
            if ((pa(n,j).gt.0.d0).and.(ma(j).le.mina)) mina = ma(j)
            if ((pb(n,j).gt.0.d0).and.(mb(j).ge.maxb)) maxb = mb(j)
            if ((pb(n,j).gt.0.d0).and.(mb(j).ge.minb)) minb = mb(j)
            if (pa(n,j).ge.pa0) then
               pa0 = pa(n,j)
               jpa = j
            endif
            if (pb(n,j).ge.pb0) then
               pb0 = pb(n,j)
               jpb = j
            endif
            if((pa(n,j).gt.0.d0).or.(pb(n,j).gt.0.d0)) jp = j
   10    continue
         ymax(n) = (maxa - maxb)/15.d0
         ymin(n) = (mina - minb)/15.d0
         ymp(n) = (ma(jpa) - mb(jpb))/15.d0
         ypm(n) = pa(n,jp) - pb(n,jp)
   50 continue
      return
      end
c***********************************************************************
      subroutine fpx(gam,u5,u10,px)
      implicit real*8 (a-h,o-z)
      parameter (ndec=100,ntype=5)
      dimension px(ndec,ntype,2)
      dimension pa(4,ndec),pb(4,ndec),ma(4),mb(4)
      dimension d1(ndec),d2(ndec),d3(ndec)
      dimension ymax(ndec),ymin(ndec),ypm(ndec),ymp(ndec)
      common /mats/pa,pb,ma,mb,d1,d2,d3
      common /ys/ymax,ymin,ypm,ymp
c
      do 20 n=1,ndec
         yeu = d1(n)*u5 + d2(n)*u10 + d3(n)
         ex1 = dexp(-gam*yeu)
         px(n,1,1) = 1.d0/(1.d0+ex1)
         ex2 = dexp(-gam*ymax(n))
         px(n,2,1) = 1.d0/(1.d0+ex2)
         ex3 = dexp(-gam*ymin(n))
         px(n,3,1) = 1.d0/(1.d0+ex3)
         ex4 = dexp(-gam*ypm(n))
         px(n,4,1) = 1.d0/(1.d0+ex4)
         ex5 = dexp(-gam*ymp(n))
         px(n,5,1) = 1.d0/(1.d0+ex5)
         do 10 k=1,ntype
            px(n,k,2) = 1.d0 - px(n,k,1)
   10    continue
   20 continue
c
      return
      end
c***********************************************************************
      double precision function gamoeb1(x)
      implicit real*8 (a-h,o-z)
      parameter (np=7,ndec=100,ntype=5)
      dimension x(np),px(ndec,ntype,2),a(ntype)
      dimension is1(ndec)
      common /ch/is1
c
      gam=1.2d2/(1.d0+dexp(x(1)))
      u5=1.d0/(1.d0+dexp(x(2)))
      u10 = u5 + (1.d0-u5)/(1.d0+dexp(x(3)))
      at = 1.d0
      do 5 k=1,ntype-1
         a(k)=at/(1.d0+dexp(x(k+3)))
         at = at - a(k)
    5 continue
      if (at.gt.0.d0) then
         a(ntype) = at
      else
         a(ntype) = 0.d0
      endif
c
      call fpx(gam,u5,u10,px)
c
      pt=1.d0
      do 20 n=1,75
         j = is1(n)
         pn = 0.d0
         do 10 k=1,ntype
            pn = pn + a(k)*px(n,k,j)
   10    continue
         pt = pt*pn
   20 continue
      gamoeb1 = -dlog(pt)
c
      return
      end
c***********************************************************************
      double precision function gamoeb2(x)
      implicit real*8 (a-h,o-z)
      parameter (np=7,ndec=100,ntype=5)
      dimension x(np),px(ndec,ntype,2),a(ntype)
      dimension is1(ndec)
      common /ch/is1
c
      gam=1.2d2/(1.d0+dexp(x(1)))
      u5=1.d0/(1.d0+dexp(x(2)))
      u10 = u5 + (1.d0-u5)/(1.d0+dexp(x(3)))
      at = 1.d0
      do 5 k=1,ntype-1
         a(k)=at/(1.d0+dexp(x(k+3)))
         at = at - a(k)
    5 continue
      if (at.gt.0.d0) then
         a(ntype) = at
      else
         a(ntype) = 0.d0
      endif
c
      call fpx(gam,u5,u10,px)
c
      pt=1.d0
      do 20 n=76,ndec
         j = is1(n)
         pn = 0.d0
         do 10 k=1,ntype
            pn = pn + a(k)*px(n,k,j)
   10    continue
         pt = pt*pn
   20 continue
      if (pt .gt. 1.d-300) then
        gamoeb2 = -dlog(pt)
      else
        gamoeb2 = 7.d2
      endif
c
      return
      end
c***********************************************************************
      subroutine jnit(xint,p,y,delta,gamoeb)
      implicit real*8 (a-h,o-z) 
      parameter (np=7,mp=np+1)
      dimension q(np,np),p(mp,np),x(np),y(mp) 
      dimension xint(np),xx(np),xxx(np)
      external gamoeb 
c 
      do 1 k=1,np
         xxx(k)=xint(k)
    1 continue
c
c      write(*,'(//''Initial Parameter Values:'')')
c      do 14 k=1,np
c         write(*,'(3x,g12.6,3x,g12.6)')xint(k),xxx(k)
c   14 continue
c
c      y0=-gamoeb(xxx)
c      write(*,'(/''initial log-likelihood = '',g14.7/)')y0
c      stop
c
      p(1,1)=1.d0
      p(2,1)=-1.d0
      if (np.eq.1) goto 51
      do 5 k=2,np
         p(1,k)=0.d0
         p(2,k)=0.d0
    5 continue
c
      do 50 it=2,np
         f0=float(it)
         f1=(dsqrt(f0*f0-1))/f0
         f2=-1.d0/f0
c
         do 11 j=1,it
            do 10 k=1,np
               q(j,k)=p(j,k)
   10       continue
   11    continue
         do 15 k=1,np
            p(1,k)=0.d0
   15    continue
         p(1,it)=1.d0
         do 21 j=2,it+1
            p(j,it)=f2
            do 20 k=1,it-1
               p(j,k)=f1*q(j-1,k)
   20       continue
   21    continue
   50 continue
   51 continue
c
c Initial step size
c      delta=0.05d0
c      write(*,'(''  stepsize = '',g12.3//)')delta
c 
      do 80 j=1,mp
         do 79 k=1,np
            p(j,k) = xxx(k) + delta*p(j,k)
   79    continue
   80 continue
c 
c Initialize y. 
      do 92 i=1,mp 
         do 91 j=1,np 
            x(j)=p(i,j) 
   91       continue 
         y(i)=gamoeb(x) 
c     write(*,'(i2,9(1x,g12.4),g14.7)')i,(x(j),j=1,np),y(i)
   92    continue 
      return 
      end 
c***********************************************************************
      subroutine genx(dseed,gam,u1,u2)
      implicit real*8 (a-h,o-z)
      parameter (ndec=100)
      dimension px(ndec,2)
      dimension is1(ndec)
      common /ch/is1
c
      call feut(gam,u1,u2,px)
c
      do 20 n=1,ndec
         pn = px(n,1)
         t1 = ggubfs(dseed)
         if (t1 .le. pn) then
            is1(n) = 1
         else
            is1(n) = 2
         endif
c         write(*,'(2i3)')is1(n)
   20 continue
      return
      end
c========================================================================
      double precision function ggubfs (dseed)
      double precision   dseed
      double precision   d2p31m,d2p31
      data              d2p31m/2147483647.d0/
      data              d2p31 /2147483648.d0/
      dseed = dmod(16807.d0*dseed,d2p31m)
      ggubfs = dseed / d2p31
      return
      end
c*********************************************************************
c ************** newlibrary of subroutines ***************************
      subroutine avg(p,x,np,mp)
      implicit real*8 (a-h,o-z)
      dimension p(mp,np),x(np)
c
c compute average values of vertices of the simplex.
      do 15 k=1,np
         do 14 j=2,mp
   14       p(1,k)=p(1,k)+p(j,k)
   15    x(k)=p(1,k)/mp
      return
      end
c------------------------------------------------------------------------
      subroutine amoeba(p,y,mp,np,ndim,ftol,funk,iter)
      implicit real*8 (a-h,o-z)
      parameter (nmax=20,alpha=1.d0,beta=0.5d0,gamma=2.d0)
      dimension p(mp,np),y(mp),pr(nmax),prr(nmax),pbar(nmax)
c
c This is the simplex algorithm from Numerical Recipes, with 
c modifications.
c When amoeba is called, iter must be set equal to the maximum
c number of iterations that will be allowed.  On return,
c iter equals the actual number of iterations.
c Thus, an example of a calling sequence would be:
c          itmax=5000
c          iter=5000
c          call aboeba(p,y,mp,np,ndim,ftol,funk,iter)
c          if (iter.ge.itmax) ...
c On return, if the if condition evaluates to true, convergence
c has not been acheived in the alloted number of iterations.
c     --PWW   1/26/94
c
c
      itmax=iter
      if (itmax.le.0) write(*,100)itmax
  100 format('*****WARNING: Maximum no. of iterations passed to',
     &       ' amoeba=',i6,'*****')
      mpts=ndim+1
      iter=0
    1 ilo=1
      if(y(1).gt.y(2))then
        ihi=1
        inhi=2
      else
        ihi=2
        inhi=1
      endif
      do 11 i=1,mpts
        if(y(i).lt.y(ilo)) ilo=i
        if(y(i).gt.y(ihi))then
        inhi=ihi
        ihi=i
        else if(y(i).gt.y(inhi))then
        if(i.ne.ihi) inhi=i
        endif
   11 continue
      rtol=2.d0*dabs(y(ihi)-y(ilo))/(dabs(y(ihi))+dabs(y(ilo)))
      if(rtol.lt.ftol)return
c
c check number of iterations.
      if (iter.ge.itmax) then
         write(*,101)itmax
         return
         end if
  101 format('*****Maximum iterations (',i6,') exceeded in amoeba.'/
     &       '     Control passed to calling routine.*****')
      iter=iter+1
      do 12 j=1,ndim
        pbar(j)=0.d0
   12 continue
      do 14 i=1,mpts
        if(i.ne.ihi)then
        do 13 j=1,ndim
        pbar(j)=pbar(j)+p(i,j)
   13   continue
        endif
   14 continue
      do 15 j=1,ndim
        pbar(j)=pbar(j)/ndim
        pr(j)=(1.d0+alpha)*pbar(j)-alpha*p(ihi,j)
   15 continue
      ypr=funk(pr)
      if(ypr.le.y(ilo))then
        do 16 j=1,ndim
        prr(j)=gamma*pr(j)+(1.d0-gamma)*pbar(j)
   16   continue
        yprr=funk(prr)
        if(yprr.lt.y(ilo))then
        do 17 j=1,ndim
        p(ihi,j)=prr(j)
   17   continue
        y(ihi)=yprr
        else
        do 18 j=1,ndim
        p(ihi,j)=pr(j)
   18   continue
        y(ihi)=ypr
        endif
      else if(ypr.ge.y(inhi))then
        if(ypr.lt.y(ihi))then
        do 19 j=1,ndim
        p(ihi,j)=pr(j)
   19   continue
        y(ihi)=ypr
        endif
        do 21 j=1,ndim
        prr(j)=beta*p(ihi,j)+(1.d0-beta)*pbar(j)
   21   continue
        yprr=funk(prr)
        if(yprr.lt.y(ihi))then
        do 22 j=1,ndim
        p(ihi,j)=prr(j)
   22   continue
        y(ihi)=yprr
        else
        do 24 i=1,mpts
        if(i.ne.ilo)then
        do 23 j=1,ndim
                pr(j)=0.5d0*(p(i,j)+p(ilo,j))
                p(i,j)=pr(j)
   23   continue
        y(i)=funk(pr)
        endif
   24   continue
        endif
      else
        do 25 j=1,ndim
        p(ihi,j)=pr(j)
   25   continue
        y(ihi)=ypr
      endif
      go to 1
      end
c======================================================================== 
      subroutine sort(n,ra)
      implicit real*8 (a-h,o-z)
      dimension ra(n)
      l=n/2+1
      ir=n
10    continue
        if(l.gt.1)then
        l=l-1
        rra=ra(l)
        else
        rra=ra(ir)
        ra(ir)=ra(1)
        ir=ir-1
        if(ir.eq.1)then
        ra(1)=rra
        return
        endif
        endif
        i=l  
        j=l+l
20      if(j.le.ir)then
        if(j.lt.ir)then
        if(ra(j).lt.ra(j+1))j=j+1
        endif
        if(rra.lt.ra(j))then
        ra(i)=ra(j)
        i=j
        j=j+j
        else
        j=ir+1
        endif
        go to 20
        endif
        ra(i)=rra
      go to 10
      end      
c***************************************************************************
      subroutine fcn1(n,x,f)
      integer n
      double precision  x(n), f, famoeb1
      external famoeb1
      f = -famoeb1(x)
      return
      end
c***************************************************************************
      subroutine gcn1(n,x,g)
      integer n
      double precision  x(n), g, gamoeb1
      external gamoeb1
      g = -gamoeb1(x)
      return
      end
c***************************************************************************
C  Simulated Annealing Version: 3.2
C  Date: 1/22/94.
C
      SUBROUTINE SA(N,X,MAX,RT,EPS,NS,NT,NEPS,MAXEVL,LB,UB,C,IPRINT,
     1              ISEED1,ISEED2,T,VM,XOPT,FOPT,NACC,NFCNEV,NOBDS,IER,
     2              FSTAR,XP,NACP,FCN)

C  Type all external variables.
      DOUBLE PRECISION  X(N), LB(N), UB(N), C(N), VM(N), FSTAR(NEPS),
     1                  XOPT(N), XP(N), T, EPS, RT, FOPT
      INTEGER  NACP(N), N, NS, NT, NEPS, NACC, MAXEVL, IPRINT,
     1         NOBDS, IER, NFCNEV, ISEED1, ISEED2
      LOGICAL  MAX

C  Type all internal variables.
      DOUBLE PRECISION  F, FP, P, PP, RATIO
      INTEGER  NUP, NDOWN, NREJ, NNEW, LNOBDS, H, I, J, M
      LOGICAL  QUIT

C  Type all functions.
      DOUBLE PRECISION  EXPREP
      REAL  RANMAR

C  Initialize the random number generator RANMAR.
      CALL RMARIN(ISEED1,ISEED2)

C  Set initial values.
      NACC = 0
      NOBDS = 0
      NFCNEV = 0
      IER = 99

      DO 10, I = 1, N
         XOPT(I) = X(I)
         NACP(I) = 0
10    CONTINUE

      DO 20, I = 1, NEPS
         FSTAR(I) = 1.0D+20
20    CONTINUE 

C  If the initial temperature is not positive, notify the user and 
C  return to the calling routine. 
      IF (T .LE. 0.0) THEN
         WRITE(*,'(/,''  THE INITIAL TEMPERATURE IS NOT POSITIVE. ''
     1             /,''  RESET THE VARIABLE T. ''/)')
         IER = 3
         RETURN
      END IF

C  If the initial value is out of bounds, notify the user and return
C  to the calling routine.
      DO 30, I = 1, N
         IF ((X(I) .GT. UB(I)) .OR. (X(I) .LT. LB(I))) THEN
            CALL PRT1
            IER = 2
            RETURN
         END IF
30    CONTINUE

C  Evaluate the function with input X and return value as F.
      CALL FCN(N,X,F)
C  If the function is to be minimized, switch the sign of the function.
C  Note that all intermediate and final output switches the sign back
C  to eliminate any possible confusion for the user.
      IF(.NOT. MAX) F = -F
      NFCNEV = NFCNEV + 1
      FOPT = F
      FSTAR(1) = F
      IF(IPRINT .GE. 1) CALL PRT2(MAX,N,X,F)

C  Start the main loop. Note that it terminates if (i) the algorithm
C  succesfully optimizes the function or (ii) there are too many
C  function evaluations (more than MAXEVL).
100   NUP = 0
      NREJ = 0
      NNEW = 0
      NDOWN = 0
      LNOBDS = 0

      DO 400, M = 1, NT
         DO 300, J = 1, NS
            DO 200, H = 1, N

C  Generate XP, the trial value of X. Note use of VM to choose XP.
               DO 110, I = 1, N
                  IF (I .EQ. H) THEN
                     XP(I) = X(I) + (RANMAR()*2.- 1.) * VM(I)
                  ELSE
                     XP(I) = X(I)
                  END IF

C  If XP is out of bounds, select a point in bounds for the trial.
                  IF((XP(I) .LT. LB(I)) .OR. (XP(I) .GT. UB(I))) THEN
                    XP(I) = LB(I) + (UB(I) - LB(I))*RANMAR()
                    LNOBDS = LNOBDS + 1
                    NOBDS = NOBDS + 1
                    IF(IPRINT .GE. 3) CALL PRT3(MAX,N,XP,X,FP,F)
                  END IF
110            CONTINUE

C  Evaluate the function with the trial point XP and return as FP.
               CALL FCN(N,XP,FP)

               IF(.NOT. MAX) FP = -FP
               NFCNEV = NFCNEV + 1
               IF(IPRINT .GE. 3) CALL PRT4(MAX,N,XP,X,FP,F)

C  If too many function evaluations occur, terminate the algorithm.
               IF(NFCNEV .GE. MAXEVL) THEN
                  CALL PRT5
                  IF (.NOT. MAX) FOPT = -FOPT
                  IER = 1
                  RETURN
               END IF

C  Accept the new point if the function value increases.
               IF(FP .GE. F) THEN
                  IF(IPRINT .GE. 3) THEN
                     WRITE(*,'(''  POINT ACCEPTED'')')
                  END IF
                  DO 120, I = 1, N
                     X(I) = XP(I)
120               CONTINUE
                  F = FP
                  NACC = NACC + 1
                  NACP(H) = NACP(H) + 1
                  NUP = NUP + 1

C  If greater than any other point, record as new optimum.
                  IF (FP .GT. FOPT) THEN
                  IF(IPRINT .GE. 3) THEN
                     WRITE(*,'(''  NEW OPTIMUM'')')
                  END IF
                     DO 130, I = 1, N
                        XOPT(I) = XP(I)
130                  CONTINUE
                     FOPT = FP
                     NNEW = NNEW + 1
                  END IF

C  If the point is lower, use the Metropolis criteria to decide on
C  acceptance or rejection.
               ELSE
                  P = EXPREP((FP - F)/T)
                  PP = RANMAR()
                  IF (PP .LT. P) THEN
                     IF(IPRINT .GE. 3) CALL PRT6(MAX)
                     DO 140, I = 1, N
                        X(I) = XP(I)
140                  CONTINUE
                     F = FP
                     NACC = NACC + 1
                     NACP(H) = NACP(H) + 1
                     NDOWN = NDOWN + 1
                  ELSE
                     NREJ = NREJ + 1
                     IF(IPRINT .GE. 3) CALL PRT7(MAX)
                  END IF
               END IF

200         CONTINUE
300      CONTINUE

C  Adjust VM so that approximately half of all evaluations are accepted.
         DO 310, I = 1, N
            RATIO = DFLOAT(NACP(I)) /DFLOAT(NS)
            IF (RATIO .GT. .6) THEN
               VM(I) = VM(I)*(1. + C(I)*(RATIO - .6)/.4)
            ELSE IF (RATIO .LT. .4) THEN
               VM(I) = VM(I)/(1. + C(I)*((.4 - RATIO)/.4))
            END IF
            IF (VM(I) .GT. 0.5d0*(UB(I)-LB(I))) THEN
               VM(I) = 0.5d0*(UB(I) - LB(I))
            END IF
310      CONTINUE

         IF(IPRINT .GE. 2) THEN
            CALL PRT8(N,VM,XOPT,X)
         END IF

         DO 320, I = 1, N
            NACP(I) = 0
320      CONTINUE

400   CONTINUE

      IF(IPRINT .GE. 1) THEN
         CALL PRT9(MAX,N,T,XOPT,VM,FOPT,NUP,NDOWN,NREJ,LNOBDS,NNEW)
      END IF

C  Check termination criteria.
      QUIT = .FALSE.
      FSTAR(1) = F
      IF ((FOPT - FSTAR(1)) .LE. EPS) QUIT = .TRUE.
      DO 410, I = 1, NEPS
         IF (ABS(F - FSTAR(I)) .GT. EPS) QUIT = .FALSE.
410   CONTINUE

C  Terminate SA if appropriate.
      IF (QUIT) THEN
         DO 420, I = 1, N
            X(I) = XOPT(I)
420      CONTINUE
         IER = 0
         IF (.NOT. MAX) FOPT = -FOPT
         IF(IPRINT .GE. 1) CALL PRT10
         RETURN
      END IF

C  If termination criteria is not met, prepare for another loop.
      NT = NT-5
      if (nt.eq.0) then
c         write(*,'("***** NT reached 0 ******")')
         if (.not. max) fopt = -fopt
         return
      endif
      T = RT*T
      DO 430, I = NEPS, 2, -1
         FSTAR(I) = FSTAR(I-1)
430   CONTINUE
      F = FOPT
      DO 440, I = 1, N
         X(I) = XOPT(I)
440   CONTINUE

C  Loop again.
      GO TO 100

      END

      FUNCTION  EXPREP(RDUM)
C  This function replaces exp to avoid under- and overflows and is
C  designed for IBM 370 type machines. It may be necessary to modify
C  it for other machines. Note that the maximum and minimum values of
C  EXPREP are such that they has no effect on the algorithm.

      DOUBLE PRECISION  RDUM, EXPREP

      IF (RDUM .GT. 174.) THEN
         EXPREP = 3.69D+75
      ELSE IF (RDUM .LT. -180.) THEN
         EXPREP = 0.0
      ELSE
         EXPREP = EXP(RDUM)
      END IF

      RETURN
      END

      subroutine RMARIN(IJ,KL)
C  This subroutine and the next function generate random numbers. See
C  the comments for SA for more information. The only changes from the
C  orginal code is that (1) the test to make sure that RMARIN runs first
C  was taken out since SA assures that this is done (this test didn't
C  compile under IBM's VS Fortran) and (2) typing ivec as integer was
C  taken out since ivec isn't used. With these exceptions, all following
C  lines are original.

C This is the initialization routine for the random number generator
C     RANMAR()
C NOTE: The seed variables can have values between:    0 <= IJ <= 31328
C                                                      0 <= KL <= 30081
      real U(97), C, CD, CM
      integer I97, J97
      common /raset1/ U, C, CD, CM, I97, J97
      if( IJ .lt. 0  .or.  IJ .gt. 31328  .or.
     *    KL .lt. 0  .or.  KL .gt. 30081 ) then
          print '(A)', ' The first random number seed must have a value
     *between 0 and 31328'
          print '(A)',' The second seed must have a value between 0 and
     *30081'
            stop
      endif
      i = mod(IJ/177, 177) + 2
      j = mod(IJ    , 177) + 2
      k = mod(KL/169, 178) + 1
      l = mod(KL,     169)
      do 2 ii = 1, 97
         s = 0.0
         t = 0.5
         do 3 jj = 1, 24
            m = mod(mod(i*j, 179)*k, 179)
            i = j
            j = k
            k = m
            l = mod(53*l+1, 169)
            if (mod(l*m, 64) .ge. 32) then
               s = s + t
            endif
            t = 0.5 * t
3        continue
         U(ii) = s
2     continue
      C = 362436.0 / 16777216.0
      CD = 7654321.0 / 16777216.0
      CM = 16777213.0 /16777216.0
      I97 = 97
      J97 = 33
      return
      end

      function ranmar()
      real U(97), C, CD, CM
      integer I97, J97
      common /raset1/ U, C, CD, CM, I97, J97
         uni = U(I97) - U(J97)
         if( uni .lt. 0.0 ) uni = uni + 1.0
         U(I97) = uni
         I97 = I97 - 1
         if(I97 .eq. 0) I97 = 97
         J97 = J97 - 1
         if(J97 .eq. 0) J97 = 97
         C = C - CD
         if( C .lt. 0.0 ) C = C + CM
         uni = uni - C
         if( uni .lt. 0.0 ) uni = uni + 1.0
         RANMAR = uni
      return
      END

      SUBROUTINE PRT1
C  This subroutine prints intermediate output, as does PRT2 through
C  PRT10. Note that if SA is minimizing the function, the sign of the
C  function value and the directions (up/down) are reversed in all
C  output to correspond with the actual function optimization. This
C  correction is because SA was written to maximize functions and
C  it minimizes by maximizing the negative a function.

      WRITE(*,'(/,''  THE STARTING VALUE (X) IS OUTSIDE THE BOUNDS ''
     1          /,''  (LB AND UB). EXECUTION TERMINATED WITHOUT ANY''
     2          /,''  OPTIMIZATION. RESPECIFY X, UB OR LB SO THAT  ''
     3          /,''  LB(I) .LT. X(I) .LT. UB(I), I = 1, N. ''/)')

      RETURN
      END

      SUBROUTINE PRT2(MAX,N,X,F)

      DOUBLE PRECISION  X(*), F
      INTEGER  N
      LOGICAL  MAX

c      WRITE(*,'(''  '')')
      CALL PRTVEC(X,N,'INITIAL X')
      IF (MAX) THEN
         WRITE(*,'(/''  INITIAL F:  '', G25.18)') F
      ELSE
         WRITE(*,'(/''  INITIAL F:  '', G25.18)') -F
      END IF

      RETURN
      END

      SUBROUTINE PRT3(MAX,N,XP,X,FP,F)

      DOUBLE PRECISION  XP(*), X(*), FP, F
      INTEGER  N
      LOGICAL  MAX

      WRITE(*,'(''  '')')
      CALL PRTVEC(X,N,'CURRENT X')
      IF (MAX) THEN
         WRITE(*,'(''  CURRENT F: '',G25.18)') F
      ELSE
         WRITE(*,'(''  CURRENT F: '',G25.18)') -F
      END IF
      CALL PRTVEC(XP,N,'TRIAL X')
      WRITE(*,'(''  POINT REJECTED SINCE OUT OF BOUNDS'')')

      RETURN
      END

      SUBROUTINE PRT4(MAX,N,XP,X,FP,F)

      DOUBLE PRECISION  XP(*), X(*), FP, F
      INTEGER  N
      LOGICAL  MAX

      WRITE(*,'(''  '')')
      CALL PRTVEC(X,N,'CURRENT X')
      IF (MAX) THEN
         WRITE(*,'(''  CURRENT F: '',G25.18)') F
         CALL PRTVEC(XP,N,'TRIAL X')
         WRITE(*,'(''  RESULTING F: '',G25.18)') FP
      ELSE
         WRITE(*,'(''  CURRENT F: '',G25.18)') -F
         CALL PRTVEC(XP,N,'TRIAL X')
         WRITE(*,'(''  RESULTING F: '',G25.18)') -FP
      END IF

      RETURN
      END

      SUBROUTINE PRT5

      WRITE(*,'(/,''  TOO MANY FUNCTION EVALUATIONS; CONSIDER ''
     1          /,''  INCREASING MAXEVL OR EPS, OR DECREASING ''
     2          /,''  NT OR RT. THESE RESULTS ARE LIKELY TO BE ''
     3          /,''  POOR.'',/)')

      RETURN
      END

      SUBROUTINE PRT6(MAX)

      LOGICAL  MAX

      IF (MAX) THEN
         WRITE(*,'(''  THOUGH LOWER, POINT ACCEPTED'')')
      ELSE
         WRITE(*,'(''  THOUGH HIGHER, POINT ACCEPTED'')')
      END IF

      RETURN
      END

      SUBROUTINE PRT7(MAX)

      LOGICAL  MAX

      IF (MAX) THEN
         WRITE(*,'(''  LOWER POINT REJECTED'')')
      ELSE
         WRITE(*,'(''  HIGHER POINT REJECTED'')')
      END IF

      RETURN
      END

      SUBROUTINE PRT8(N,VM,XOPT,X)

      DOUBLE PRECISION  VM(*), XOPT(*), X(*)
      INTEGER  N

      WRITE(*,'(/,
     1  '' INTERMEDIATE RESULTS AFTER STEP LENGTH ADJUSTMENT'',/)')
      CALL PRTVEC(VM,N,'NEW STEP LENGTH (VM)')
      CALL PRTVEC(XOPT,N,'CURRENT OPTIMAL X')
      CALL PRTVEC(X,N,'CURRENT X')
      WRITE(*,'('' '')')

      RETURN
      END

      SUBROUTINE PRT9(MAX,N,T,XOPT,VM,FOPT,NUP,NDOWN,NREJ,LNOBDS,NNEW)

      DOUBLE PRECISION  XOPT(*), VM(*), T, FOPT
      INTEGER  N, NUP, NDOWN, NREJ, LNOBDS, NNEW, TOTMOV
      LOGICAL  MAX

      TOTMOV = NUP + NDOWN + NREJ

      WRITE(*,'(/,
     1  '' INTERMEDIATE RESULTS BEFORE NEXT TEMPERATURE REDUCTION'',/)')
      WRITE(*,'(''  CURRENT TEMPERATURE:            '',G12.5)') T
      IF (MAX) THEN
         WRITE(*,'(''  MAX FUNCTION VALUE SO FAR:  '',G25.18)') FOPT
         WRITE(*,'(''  TOTAL MOVES:                '',I8)') TOTMOV
         WRITE(*,'(''     UPHILL:                  '',I8)') NUP
         WRITE(*,'(''     ACCEPTED DOWNHILL:       '',I8)') NDOWN
         WRITE(*,'(''     REJECTED DOWNHILL:       '',I8)') NREJ
         WRITE(*,'(''  OUT OF BOUNDS TRIALS:       '',I8)') LNOBDS
         WRITE(*,'(''  NEW MAXIMA THIS TEMPERATURE:'',I8)') NNEW
      ELSE
         WRITE(*,'(''  MIN FUNCTION VALUE SO FAR:  '',G25.18)') -FOPT
         WRITE(*,'(''  TOTAL MOVES:                '',I8)') TOTMOV
         WRITE(*,'(''     DOWNHILL:                '',I8)')  NUP
         WRITE(*,'(''     ACCEPTED UPHILL:         '',I8)')  NDOWN
         WRITE(*,'(''     REJECTED UPHILL:         '',I8)')  NREJ
         WRITE(*,'(''  TRIALS OUT OF BOUNDS:       '',I8)')  LNOBDS
         WRITE(*,'(''  NEW MINIMA THIS TEMPERATURE:'',I8)')  NNEW
      END IF
      CALL PRTVEC(XOPT,N,'CURRENT OPTIMAL X')
      CALL PRTVEC(VM,N,'STEP LENGTH (VM)')
      WRITE(*,'('' '')')

      RETURN
      END

      SUBROUTINE PRT10

      WRITE(*,'(/,''  SA ACHIEVED TERMINATION CRITERIA. IER = 0. '',
     &/)')

      RETURN
      END

      SUBROUTINE PRTVEC(VECTOR,NCOLS,NAME)
C  This subroutine prints the double precision vector named VECTOR.
C  Elements 1 thru NCOLS will be printed. NAME is a character variable
C  that describes VECTOR. Note that if NAME is given in the call to
C  PRTVEC, it must be enclosed in quotes. If there are more than 10
C  elements in VECTOR, 10 elements will be printed on each line.

      INTEGER NCOLS
      DOUBLE PRECISION VECTOR(NCOLS)
      CHARACTER *(*) NAME

      WRITE(*,1001) NAME

      IF (NCOLS .GT. 6) THEN
         LINES = INT(NCOLS/6.)

         DO 100, I = 1, LINES
            LL = 6*(I - 1)
            WRITE(*,1000) (VECTOR(J),J = 1+LL, 6+LL)
  100    CONTINUE

         WRITE(*,1000) (VECTOR(J),J = 7+LL, NCOLS)
      ELSE
         WRITE(*,1000) (VECTOR(J),J = 1, NCOLS)
      END IF

 1000 FORMAT( 10(G12.5,1X))
 1001 FORMAT(/,25X,A)

      RETURN
      END
